Traveling waves of infection in the Hantavirus epidemics 
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Traveling waves are analyzed in a model of the Hantavirus infection of deer mice. The existence 
of two kinds of wave phenomena is predicted. An environmental parameter governs a transition 
between two regimes of propagation. In one of them the front of infection lags behind at a constant 
rate. In the other, fronts of susceptible and infected mice travel at the same speed, separated by 
a constant delay. The dependence of the delay on system parameters is analyzed numerically and 
through a piecewise linearization. 

PACS numbers: 87.19.Xx, 87.23.Cc, 05.45.-a 



I. INTRODUCTION 

The Hantavirus zoonosis is a persistent problem in 
many regions of the world 0. Each species of Han- 
tavirus is almost exclusively associated with a single ro- 
dent species from which, eventually, humans may result 
infected. The disease can range from mild to very se- 
vere, with a mortality in excess of 50%. Such is the case 
of the Hantavirus Pulmonary Syndrome (HPS) produced 
by the virus associated with the deer mouse {Peromyscus 
maniculatus) . In the North American Southwest, indeed, 
a serious outbreak of this disease in 1993 led to the iden- 
tification of the virus and its association with the deer 
mouse 1^ . Since then, enormous effort and resources have 
been devoted to the understanding of the ecology and the 
epidemiology of the virus-mouse association, with the ul- 
timate goal of being able to predict and prevent the risk 
for the human population ^ . 

As observed by Yates et al. [D , the role of the envi- 
ronment seems to be determinant in the prevalence of 
the infection within the mouse population, and the re- 
lated incidence of the HPS. Its has been observed that 
the disease can disappear completely from a local popu- 
lation during times of adverse environmental conditions, 
only to reappear sporadically ^. Besides, there are 
indications of focality of the infection in "refugia" of the 
mouse population. Both phenomena are most probably 
related, the refugia acting as reservoirs of the virus during 
times when the infection has disappeared from most of 
the landscape. When environmental conditions change, 
the infection spreads again from these refugia. 

In a recent work |Q , Abramson and Kenkre have shown 
that a simple epidemic model is able to display qualita- 
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tively these behaviors, as the result of a bifurcation of 
the equilibrium states of the system as controlled by the 
carrying capacity of the medium. The purpose of the 
present paper is to analyze the dynamics of simple trav- 
eling waves, as a model of the mechanisms in which an 
epidemic wave of Hantavirus might propagate into a pre- 
viously uninfected region. 



II. SPATIALLY EXTENDED MODEL 

The model of Ref . § is a mean-field continuous model 
which has been intentionally kept simple enough to facil- 
itate the comprehension of the basic mechanisms, and at 
the same time to incorporate as many known facts about 
the ecology and epidemiology of the biological system 
as possible. The reader is referred to for a detailed 
discussion, which we summarize here. It is known that 
the virus does not affect properties such as the mortality 
of the mice, so that no appreciable difference in death 
rate, for example, is to be expected between susceptible 
and infected mice. It is also not transmitted to new- 
borns, so that all mice are born susceptible. The infection 
is transmitted from mouse to mouse through individual 
contacts, presumably during fights. More general facts 
of the ecology of Peromyscus indicate that adults occa- 
sionally shift their home range to nearby locations, in 
particular if these are vacant ^ . This enables us to 
model the transport of mice as a diffusion process. Fi- 
nally, intra-species competition for resources indicate a 
saturated population growth, which has been observed 
to be of a logistic form in the laboratory ||l^. Logistic 
growth is also a well established metaphor of the dynam- 
ics of a self-limitating population pd| |. 

For the sake of simplicity, assume further that the only 
population structure is the division of the whole popula- 
tion into susceptible and infected mice, denoted by Ms 
and Mj respectively. With these ingredients, the model 
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is described by the following equations: 



dMs 
dt 

dMi 



bM - cMs 
-cMi 



MsM 
MiM 



- aMgMj + D\J^Mil) 



K{x,t) 



aMsMi + DV^Ms- (2) 



Observe that the carrying capacity K{x,t), containing 
the most direct effect of the environment on the mouse 
population, is allowed, in our model, a spatial and a tem- 
poral variation, to accommodate for a diversity of habi- 
tats and temporal phenomena. The latter comprise the 
yearly variation due to seasonality, but also sporadic fluc- 
tuations such as droughts and El Nino effects. 

The sum of the two equations reduces to a single 

equation for the whole population: 



dM 
~dt 



[h - c)M 1 - 



M 



{b-c)K 



(3) 



This is Fisher's equation, originally proposed as a de- 
terministic model of the spread of a favored gene in a 
population , and which eventually became a standard 
model for a self regulated field in a diversity of situations 

In Ref. m it was shown that, as a function of K, the 
system undergoes a bifurcation between a stable state 
with only susceptible mice (and Mj = 0) to a stable 
state with both positive populations. The value of the 
critical carrying capacity is a function of the parameters 
in the following way: 



a(b — c) 



(4) 



This critical value does not depend on D, and the same 
bifurcation is observed either in a space- independent sys- 
tem {D = 0) or in a homogeneous extended one in the 
presence of diffusion. In an inhomogeneous situation, for 
moderate values of the diffusion coefhcient, the infected 
subpopulation remains restricted to those places where 
K > Kc, becoming extinct in the rest. During times of 
adverse environmental conditions, these regions become 
isolated and constitute the observed refugia. Figure |l| 
shows a typical situation of this phenomenon. A simu- 
lated, albeit realistic, landscape of K{x) has been derived 
from satellite images in Northern Patagonia. The carry- 
ing capacity is supposed proportional to the vegetation 
cover, and results highest along a river, clearly inferred 
from the density plots. These show the distribution of 
the populations of susceptible and infected mice. It can 
be seen that susceptible mice cover most of the range. 
Meanwhile, the infected population has become extinct 
in most of it, and persists only in some of the places of 
highest K. The distributions shown in Fig. |l| constitute a 
stable equilibrium of the system, found by numerical res- 
olution of Eqs. (0)-(H) from an arbitrary initial condition, 
and with zero-current boundary conditions. 
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FIG. 1: Density plots showing characteristic distribution of 
susceptible and infected mice in an inhomogeneous landscape, 
where the carrying capacity has been modeled according to 
the vegetation, derived from satellite imagery. 



III. TRAVELING WAVES 

When conditions in the landscape change, how do the 
infected phase evolve from the refugia, retracting from 
or, more importantly, invading previously uninfected re- 
gions? this is the primary question we address in the 
present paper. Fisher's equation (|^) has found wide ap- 
plicability for the description of the formation and prop- 
agation of spatio-temporal patterns. Traveling wave so- 
lutions of Fisher's equation paradigmatically show how 
a stable phase invades an unstable one in the form of a 
front propagating at a certain speed. 

There is no reason to suppose a priori that the two 
waves, susceptible and infective, will travel at the same 
speed. Accordingly, we use an ansatz which incorporates 
two independent traveling waves. In one dimension, zi — 
X — vst in Eq. (|I]) and Z2 = x — vjt in Eq. (^. This gives 
the following second-order system of ordinary differential 
equations: 



D 



(PMsizi) , dMsizi) 



dz\ 



vs- 



dz-i 



f{Ms,Mi) = 0, (5) 
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D 



VI- 



dz2 



+ g{Ms,Mi) = 0, (6) 



where vs and vj are the speeds of the susceptible and 
infected waves respectively, and / and g are the non- 
diffusive term in (|^)-(||). 

There are two interesting scenarios for these waves. In 
the first one, a large part of the system is initially at 
a state of low carrying capacity, below Kc, and conse- 
quently the population consists of uninfected mice only, 
at the stable equilibrium. Let us suppose that this region 
is in contact with a refugium. If environmental changes 
occur, and the whole region finds itself at a value of the 
carrying capacity K > the population will be out 
of equilibrium. Two processes will occur simultaneously: 
the population of susceptible mice will evolve towards a 
new equilibrium, and a wave of infected mice will advance 
from the refugium, invading the susceptible population. 
The speed of this wave can be calculated from the sta- 
bility analysis of the equilibrium states, requiring that 
the unstable infected mice density does not oscillate be- 
low zero. This unstable equilibrium is Mg = K(b — c), 
Mj = 0, and a linear stability analysis of the system 
®"® provides the following four eigenvalues: 



Ai,: 



A3, 4 — 



-v± y^V^ + AD{h ~ c) 
2D ' 
-V ± yjv"^ + 4D[6 - aK{h ~ c)] 
2D ' 



(7) 



(8) 



The requirement that Mi{z) does not oscillate below 
imposes a restriction to the radical in Eq. (||), from which 
we find the following expression for the speed of the trav- 
eling wave: 



V > 2^D \-h^aK(h-c)\. (9) 

An example of such a wave is shown in Fig. ^ja, found by a 
numerical integration of Eqs. ([l]) and (|^) in 1 dimension. 

The second interesting scenario corresponds to a sys- 
tem which is initially empty of mice, both susceptible and 
infected. This situation is always unstable within the va- 
lidity of our simple model, but it is certainly a biological 
possibility. Consider a system with K > Kc, and with 
Ms = Mj = in almost all of its range, but in contact 
with a refugium in equilibrium. A wave of both mice 
populations will develop invading the empty region. In 
fact, the total population will display just the behavior 
of a traveling wave of Fisher's equation. This wave will 
be composed of two fronts, susceptible and infected re- 
spectively, with a delay of the latter with respect to the 
former. As before, we can find the speeds of these two 
waves from a stability analysis around the correspond- 
ing equilibria. The leading wave propagates into the null 
equilibrium, Mg = Mj = 0, to which the following eigen- 
values correspond: 



Ail, 2 



-VI ± y/vf+iDc 
2D ' 



(10) 



-vs±^vi-iD{b-c) 
"''^ = 2D • (^^^ 
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FIG. 2: Traveling waves in the one dimensional model, a) 
Infection wave invading a noninfected population, b) Nonin- 
fected mice invading an empty region, followed by an infection 
wave. Model parameters: a = 0.1, & = 1, c = 0.5, D = 1, 
A" = 30. Both waves move at the same speed v — \/2 for this 
choice of parameters. 



In this situations, we require that Ms{z) does not oscil- 



late below 0, and Eq. (Ill 



) provides the restriction on the 



speed of the susceptible front: 



vs > 2^D{h~c), 



(12) 



which is, naturally, the same result as for Fisher's equa- 
tion. The second front, developed when part of the quasi- 
stable population of susceptible mice is converted into 
infected, evolves from the equilibrium Mg — K(b — c), 
Mj — 0, as in the previous scenario. Consequently, the 
same linear stability analysis apply, and from Eq. (^ we 
find a speed analogous to Eq. (||) . The front of infected 
move behind the susceptible one at a speed: 



VI 



> 2y/ D[-b + aK{b - c)], 



(13) 



which, unlike wg, does depend on the contagion rate a 
and the carrying capacity K. Figure ^jb shows such a 
situation. The density of susceptible mice rises from zero 
and lingers near the positive unstable equilibrium before 
tending to the stable one. It is remarkable that a delay 
exists between the two fronts, even when no such effect 
was explicitly considered in the dynamics (such as incu- 
bation time, or age dependence). Such delays have been 
observed in some populations and rationalized in differ- 
ent ways (see or [B| for a synthesis). 

Even though Eqs. ("), ( |l^ and (|l^ ) give only a lower 
bound to the speed of propagation of the fronts, and 
allow a continuous of speeds above this, in real situations 
only the lower bound is expected to be observed as a 
stationary solution. Higher speeds may, however, play 
a role in transient situations whose relevance in far from 
equilibrium systems such as real mice in the wild, subject 
to a fluctuating environment, cannot be underestimated. 
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FIG. 3: Delay of the infected front, as a function of time, 
following an initial condition as described in the text. Two 
regimes are shown: Kc < K < Kq and Kq < K, separated 
by the critical case K = A'o, that behaves asymptotically as 
t(l/4). 



The different functional dependence of vs and vj on 
the parameters of the system (Eqs. (p^)-(p^)) indicates 
that two regimes are possible. Indeed, when vj < vs the 
front of infection lags behind the front of susceptible at a 
delay A that increases linearly in time: A{t) ~ (vs — vi)t. 
Elementary manipulation of Eqs. (^2|) and ( |l3| ) shows 
that this occurs whenever the carrying capacity satisfies: 



Kr<K < K, 



2b- 







a[b — c) 



(14) 



where Kq is a new critical carrying capacity. At if = Kq 
the delay becomes effectively constant. For values of 
K greater than Kq, the velocities vj and vs, calculated 
from linear considerations around the equilibria, satisfy 
vi > Vs- This regime is clearly unphysical in a station- 
ary situation, since the front of susceptible necessarily 
precedes the infected one. It could be realizable and rel- 
evant in transient situations, that will be analyzed else- 
where. From numerical resolution of the system, we can 
observe that vj vs and the delay tends to a constant 
value, whenever K > Kq. Figure ^ shows the temporal 
evolution of the delay in the two regimes K < Kq and 
K > Kq, as well as in the critical case K — Kq, where 
it is seen to increase as A ^ i^/'*. It can be seen that 
there is a transient time in both regimes, that gets pro- 
gressively longer as K approaches Kq either from above 
or from below. 

With the assumption that vs = vj = v, it is possible 
to perform a piecewise linearization of Eqs. (|^)-(||) and 
find an approximate analytical expression for the front 
shapes shown in Fig. ^ and, consequently, for the delay 
A in the stationary state. The details of the calculation 
can be found in the Appendix. The main result is the 
following expression for the delay, for the case of equal 
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FIG. 4: Delay of the infected front, as a function of K. 



speeds of the two classes of mice: 



i^(h-c)a{K-KQ) 



\og{wi_W2), (15) 



where w\ and wi are complex numbers of unit modulus 
that depend on a, 6 and c, so that the logarithm is ef- 
fectively twice the phase difference between them. When 
K —t Kq , the arguments of wi and W2 tend to tt and 
respectively, so that the following critical behavior is 
predicted by the linear approximation: 



A - VD[{b - c)a{K ~ Ko)]~^/^, when K 



Kt. 



(16) 



We have analyzed the dependence of A on if by means 
of numerical resolution of the full system. In Fig. ^ we 
plot the asymptotic value of A as a function oi K — Kq 
for a variety of system parameters. The behavior is found 
to be 



A 



D[{b - c)a{K - Kq)]- 



(17) 



for values of K immediately above Kq. The exponent a 
is parameter-independent, and its value, calculated from 
the simulations is, approximately 0.388. There is a dis- 
crepancy in the exponent found numerically in the fully 
nonlinear system and the one found in the linearized ap- 
proximation, which shows the limitation of the linearized 
solution. 

This critical behavior of the delay persists while the 
wavefronts are of the kind shown in Fig. ^p. However, if 
K continues to grow the equilibrium value of the infected 
population turns greater than the equilibrium value of 
the susceptible. There is a gradual crossover to a sit- 
uation where most of the population becomes infected. 
The value of the carrying capacity at which this happens 



can be estimated from the known equilibria, Mg 
giving: 

2b 



M 



Ki 



a(b — c) 



(18) 
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Indeed, the numerical results show that the power law 
decay of the delay as a function of K starts to flatten for 
values of i^T > Ki and reaches a minimum which shows 
the same dependence that Ki does. Greater values 
of K are probably unrealistic in the arid and semi-arid 
habitats of the Southwest. It is nevertheless interesting 
to point out that the interaction of the two fronts in this 
regime results in an increase of the delay as a function 
of K . This increase is also algebraic: 

^r^{K-Kof, (19) 

with [3 K, 0.16, as found in the numerical calculations. 



IV. CONCLUSION 

We have analyzed the propagation of traveling fronts 
in 1 dimension in a simple model of the ecology and epi- 
demiology of the Hantavirus in deer mouse. We have 
found that, when a mouse- free region is in contact with 
an infected region in equilibrium, two waves propagate 
into the empty region. The first one is a wave of sus- 
ceptible mice. A wave of infected mice propagates be- 
hind it with a certain delay. Two regimes of propagation 
exist, controlled by the environmental parameter K . If 
Kc < K < Kq, the lag between the two fronts increase 
linearly in time. Conversely, ii K > Kq, the two fronts 
propagate at the same speed and the delay depends crit- 
ically on the difference K — Kq. 

The occurrence of this double regime may be of rele- 
vance for the control of the propagation of an epidemic 
wave. Indeed, the controlled reduction of K ahead of a 
propagating wave seems the most effective mean of stop- 
ping or reducing its advance. Ideally, the carrying capac- 
ity should be reduced below Kq, to ensure the complete 
extinction of the infection. However, if such a reduction 
is not feasible, the fact that Kq > Kc provides an al- 
ternative: a reduction of the carrying capacity below Kq 
would make the wave of infection start to lag more and 
more behind the wave of healthy mice. Possible imple- 
mentations of these strategies, based on the propagation 
of waves in the presence of "barriers," will be analyzed 
in detail elsewhere. 

The existence of dynamical phenomena such as these 
traveling fronts also opens the interesting possibility of 
subjecting our predictions to experimental verification. 
Controlled experiments of front propagation could be 
possible in the Sevilleta LTER facility, that the Univer- 
sity of New Mexico operates near Socorro, NM Mea- 
surements of uncontrolled mice populations along lines 
radiating from the refugia of infection will also provide 
evidence of the propagation mechanisms. The observa- 
tion of these in real mice population will provide a valu- 
able source of data to assign realistic values to the pa- 
rameters of the mathematical model. 
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V. APPENDIX: LINEARIZED SOLUTIONS 

For the purpose of finding approximate solutions of the 
waves profiles, it is better to replace the system (^)-(|^) 
with one involving A'l and Mi instead: 



DM" + vM' + {b- c)M 



~K 



= 0, 



dm;' + vM[ + q{M)Mi - aMf = 0, 



(20) 
(21) 



where q(M{z)) = — c+aM{z) — M{z)/ K and both speeds 
are assumed equal, as observed in numerical resolutions 
for the regime K > Kq. Primes denote differentiation 
with respect to z. The reason for using this system is 
that the equation for M(z), Eq. (20), being closed, can 
be solved independently of Mj. Its solution can be used 
then into Eq. ( pl| ) as a z-dependent parameter and solve 
for Mi{z). 

Linearized solutions for the traveling waves of Fisher's 
equation ( pO| ) arc well known, and essentially consist 
of two exponentials, matched smoothly at z = 0, rep- 
resenting a front that travels to the right at a speed 
V > 2y^D{b — c). Such a function needs a further sim- 
plification in order to solve Eq. (pl|). We approximate 
M{z) with a step function discontinuous at z = 0: 
M{z) = M* ^ K{b- c) ii z <0 and M{z) = if z > 0. 
Consequently, we have that q{z) = aMj — —b+aK{b—c) 
if z < and q{z) = — c if z > 0. We also linearize on both 
sides of the center value M/(— A) — Mj /2, which defines 
A as the delay between the two waves. Figure ^ shows 
the geometry of the procedure. The linearized equation 
for All breaks into three regimes: 

DM'/ + vM'j- a{Mi - Ml) = if z < -A, (22) 
DM'/ + vM'j + aMjMi = if z e [-A, 0],(23) 
DM/ + vM'j - CMi = if < z. (24) 

The solutions of these are respectively: 
b 



Miiz) 



+ K{b~c) + aie^^ iiz <-A, (25) 



6 



M2{z) = 616^+^+626^-" 

M3(z) = cie"^^ 



where: 



D 



if z € hA,0],(26) 
if < z, (27) 



(28) 



M± = 



Z? V 



:a(if-ifo), (29) 
(30) 



and the terms that diverge in infinity have already been 
cancelled. It is clear that A > and z/ < 0, as required for 
the solution not to diverge in =poo. The real part of /i± 
is always positive, and there is an imaginary part when 

The solutions Mi, M2 and need to be matched 
smoothly at z = — A and z — 0, in the following way 
(r = b — c from now on): 



Mi(-A) 


b kr 

"2^ ^ Y' 


(31) 


M2(-A) 


b kr 

"2^ ^ Y' 


(32) 


M{(-A) 


- M^(-A), 


(33) 


M2(0) 


= M3(0), 


(34) 


M^(0) 


= M^(0). 


(35) 



These are five equations with five unknowns that can 
be solved without difficulty: 



and: 



ai : 
hi : 

62 

ci - 
A = 



-b/a + kr 

~ 6 1 



2 

(A + A^-) 

2(/i+ - 

(A + /i+) 
2(^+ - 
61 + 62, 



ikr-b/a)e''+^, 
\kr-b/a)e''-^, 



■log 



(a + m+)(m- -i^) 



(A + /i_)(/i+ -z^) 



(36) 
(37) 

(38) 
(39) 

(40) 



Expressions (|25|)-(g7|), (g8|)-(|30|) and (|36|)-(^) define 
the linearized solution for Mj. An example of such a 
front wave is shown in Fig. ^. 



Using expressions (|2^)-(p0[) in ([40[), it easy to see find: 
A^ ^, ^ _ \og{wiW2), (41) 



where 



i^{b-c)a{K-K^) 



= aT7Y = ^ ' 



(42) 
(43) 




FIG. 5: Linearized solutions. System parameters are: a — 0.3, 
6 = 1, c = 0.5, K = 1Q, 0.1. 



where (pi and 02 are respectively the arguments of the 
complex numbers A + /i_|_ and — v. It is easy to see that 
(t>i and (/>2 do not depend on I?, so that the logarithm only 
corrects the dependence of A on (b — c)a{K — Kq) . It can 
be also easily observed that wi lies in the second quad- 
rant, and W2 in the fourth, and that Ivciwi — — Imw2. 
When K —t Kq , the imaginary parts tend to zero and 
the phase difference between both tend to tt. Conse- 
quently the leading behavior of A results: 



D[ib-c)aiK-Ko)]-'/'. 



(44) 



This result from the piecewise linearization does not 
agree with that found numerically in the nonlinear sys- 
tem. 
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